<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/1999/REC-html401-19991224/loose.dtd">
<html lang="en">
<head>
	<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
	<title>Homogeneous Propagation Medium Example (k-Wave)</title>
	<link rel="stylesheet" href="kwavehelpstyle.css" type="text/css">
	<meta name="description" content="Homogeneous Propagation Medium Example.">
</head>

<body><div class="content">

<h1>Homogeneous Propagation Medium Example</h1>

<p>This example provides a simple demonstration of using k-Wave for the simulation and detection of the pressure field generated by an  initial pressure distribution within a two-dimensional homogeneous propagation medium.</p>

<p>
    <ul>
        <li><a href="matlab:edit([getkWavePath('examples') 'example_ivp_homogeneous_medium.m']);" target="_top">Open the file in the MATLAB Editor</a></li>
        <li><a href="matlab:run([getkWavePath('examples') 'example_ivp_homogeneous_medium']);" target="_top">Run the file in MATLAB</a></li>
    </ul>
</p>

<h2>Contents</h2>
<div>
	<ul>
		<li><a href="#heading2">Creating the k-space grid</a></li>
		<li><a href="#heading3">Avoiding the perfectly matched layer</a></li>
		<li><a href="#heading4">Defining the medium properties</a></li>
		<li><a href="#heading5">Defining the initial pressure distribution</a></li>
		<li><a href="#heading6">Defining the sensor mask</a></li>
		<li><a href="#heading7">Running the simulation</a></li>
	</ul>
</div>

<a name="heading2"></a>
<h2>Creating the k-space grid</h2>

<p>The simulation functions in k-Wave require four input structures. These define the properties of the computational grid, the material properties of the medium, the properties and locations of any acoustic sources, and the properties and locations of the sensor points used to record the evolution of the pressure and velocity fields over time. Starting with the computational grid, the simulations are performed on a regular Cartesian mesh (for users with a background in finite-element methods, this is analogous to a structured mesh containing identical rectangular elements). The medium discretisation is performed using <code><a href="kWaveGrid.html">kWaveGrid</a></code>. Both the total number of grid points (Nx, Ny) as well as the spacing between the grid points (dx, dy) in each Cartesian direction are used to compute the discretisation, and an object of the <code><a href="kWaveGrid.html">kWaveGrid</a></code> class is returned. As the numerical techniques used in k-Wave are based heavily on the fast Fourier transform (FFT), the simulations will be fastest when the number of grid points in each direction is given by a power of two or has small prime factors (see <code><a href="matlab: doc fft2">fft2</a></code> and <code><a href="checkFactors.html">checkFactors</a></code>).</p>

<pre class="codeinput">
<span class="comment">% create the computational grid</span>
Nx = 128;           <span class="comment">% number of grid points in the x (row) direction</span>
Ny = 128;           <span class="comment">% number of grid points in the y (column) direction</span>
dx = 0.1e-3;        <span class="comment">% grid point spacing in the x direction [m]</span>
dy = 0.1e-3;        <span class="comment">% grid point spacing in the y direction [m]</span>
kgrid = kWaveGrid(Nx, dx, Ny, dy);
</pre>

<p>The object <code>kgrid</code> contains numerous properties, including matrices of the computational wavenumbers (<code>kgrid.kx</code>, <code>kgrid.ky</code>, and <code>kgrid.k</code>) and Cartesian coordinates of the grid points (<code>kgrid.x</code>, <code>kgrid.y</code>). These properties are used by many k-Wave functions.</p>

<a name="heading3"></a>
<h2>Avoiding the perfectly matched layer</h2>

<p>When the acoustic waves reach the edge of the computational domain, they are absorbed by a special type of anisotropic absorbing boundary layer known as a perfectly matched layer (PML). The effects of the layer can be seen by watching what happens to the propagating waves as they get close to the edge of the computational domain. By default, this layer occupies a strip of 20 grid points (10 grid points in 3D) around the edge of the domain <em>inside</em> the computational domain specified using <code><a href="kWaveGrid.html">kWaveGrid</a></code>. To avoid strange effects, care must be taken not to place the source or sensor points inside this layer. Alternatively, the perfectly matched layer can be set to be <em>outside</em> the computational domain set by the user. See <a href="example_na_controlling_the_pml.html">Controlling The Absorbing Boundary Layer Example</a> for more detailed instructions on how to modify the properties and position of the perfectly matched layer.</p>

<a name="heading4"></a>
<h2>Defining the medium properties</h2>

<p>For a homogeneous medium, the sound speed is set as a scalar value in SI units (i.e., metres per second). Power law acoustic absorption can also be optionally set by assigning values to <code>medium.alpha_coeff</code> and <code>medium.alpha_power</code>. These respectively correspond to the power law coefficient or prefactor in units of dB/(MHz^y cm) and the power law exponent or power y. It is useful to note that for a homogeneous medium, the computation of the acoustic pressure is not dependent on the ambient density. Consequently, if no velocity inputs or outputs are required, <code>medium.density</code> does not need to be specified. For modelling nonlinear effects, <code>medium.BonA</code> should also be specified. In this case, k-Wave includes a number of additional nonlinear terms in the discrete governing equations.</p>

<pre class="codeinput">
<span class="comment">% define the properties of the propagation medium</span>
medium.sound_speed = 1500;  <span class="comment">% [m/s]</span>
medium.alpha_coeff = 0.75;  <span class="comment">% [dB/(MHz^y cm)]</span>
medium.alpha_power = 1.5;
</pre>

<p>The grid spacing (<code>dx</code> and <code>dy</code>) and the sound speed govern the maximum frequency that the simulation grid is able to propagate. This frequency is printed to the command line when the simulation runs, and can be calculated (for each Cartesian direction) by <code>f_max_x = medium.sound_speed/(2*dx)</code>.</p>

<a name="heading5"></a>
<h2>Defining the initial pressure distribution</h2>

<p>The initial pressure distribution is set as a matrix which contains the initial pressure values for each grid point. This matrix must be the same size as the medium discretisation defined by the computational grid (i.e., it must have <code>Nx</code> rows and <code>Ny</code> columns). Several functions are included in the toolbox for the creation of simple geometric shapes. In this example, the function <code><a href="makeDisc.html">makeDisc</a></code> is used to create an initial pressure distribution in the shape of two small discs with different diameters. This distribution is assigned to the <code>p0</code> field of the <code>source</code> structure. There are no restrictions on <code>source.p0</code> except that it must be the same size as the computational grid.</p>

<pre class="codeinput">
<span class="comment">% create initial pressure distribution using makeDisc</span>
disc_magnitude = 5; <span class="comment">% [Pa]</span>
disc_x_pos = 50;    <span class="comment">% [grid points]</span>
disc_y_pos = 50;    <span class="comment">% [grid points]</span>
disc_radius = 8;    <span class="comment">% [grid points]</span>
disc_1 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; <span class="comment">% [Pa]</span>
disc_x_pos = 80;    <span class="comment">% [grid points]</span>
disc_y_pos = 60;    <span class="comment">% [grid points]</span>
disc_radius = 5;    <span class="comment">% [grid points]</span>
disc_2 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source.p0 = disc_1 + disc_2;
</pre>

<a name="heading6"></a>
<h2>Defining the sensor mask</h2>

<p>The sensor mask defines the positions within the computational domain where the pressure field is recorded at each time-step. The sensor mask can be given in one of three ways:</p>
<ol>
<li>A binary matrix which specifies the grid points that record the data (see the <a href="example_ivp_binary_sensor_mask.html">Using A Binary Sensor Mask Example</a>).</li>
<li>The grid coordinates of two opposing corners of a line (1D), rectangle (2D) or cuboid (3D) (see the <a href="example_ivp_opposing_corners_sensor_mask.html">Defining A Sensor Mask By Opposing Corners Example</a>).</li>
<li>A set of Cartesian coordinates lying within the dimensions of the computational domain.</li>
</ol>
<p>In 2D, Cartesian coordinates must be specified as a 2 x N matrix, where the Cartesian origin is assumed to be in the center of the grid. Here the function <code><a href="makeCartCircle.html">makeCartCircle</a></code> is used to set the sensor mask to the Cartesian coordinates of a set of evenly spaced points on a circle.</p>

<pre class="codeinput">
<span class="comment">% define a centered circular sensor</span>
sensor_radius = 4e-3;   <span class="comment">% [m]</span>
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);
</pre>

<p>A plot of the initial pressure distribution and the sensor mask is given below. By default, the pressure at the Cartesian points in 2D is computed at each time step using linear interpolation.</p>

<img vspace="5" hspace="5" src="images/example_ivp_homogeneous_medium_01.png" style="width:560px;height:420px;" alt="">

<a name="heading7"></a>
<h2>Running the simulation</h2>

<p>The computation is started by calling the function <code><a href="kspaceFirstOrder2D.html">kspaceFirstOrder2D</a></code> with the four input structures defined above. By default, a visualisation of the propagating wave-field and a status bar are displayed, with frame updates every ten time-steps. The default k-Wave color map displays zero as white, positive pressures as yellows to reds to black, and negative pressures as light to dark blue-greys (see <code><a href="getColorMap.html">getColorMap</a></code>). 

<pre class="codeinput">
<span class="comment">% run the simulation</span>
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);
</pre>

<p>As the function runs, status updates and computational parameters are printed to the command line.</p>

<pre class="codeinput">
Running k-Wave simulation...
  start time: 22-Feb-2017 13:33:04
  reference sound speed: 1500m/s
  dt: 20ns, t_end: 12.06us, time steps: 604
  input grid size: 128 by 128 grid points (12.8 by 12.8mm)
  maximum supported frequency: 7.5MHz
  smoothing p0 distribution...
  calculating Delaunay triangulation...
  precomputation completed in 0.67051s
  starting time loop...
  estimated simulation time 3.2012s...
  simulation completed in 3.5675s
  total computation time 4.305s
</pre>

<p>When the time loop has completed, the function returns the recorded time series at each of the sensor points defined by <code>sensor.mask</code>. This is indexed as <code>sensor_data(sensor_point_index, time_index)</code>. For a Cartesian sensor mask, the time series are returned in the same order as the Cartesian coordinates specified in <code>sensor.mask</code>. A visualisation of the sensor data recorded in this example is given below.</p>

<pre class="codeinput">
<span class="comment">% plot the simulated sensor data</span>
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;
</pre>

<img vspace="5" hspace="5" src="images/example_ivp_homogeneous_medium_02.png" style="width:560px;height:420px;" alt="">

</div></body></html>